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Abstract 

Following the usual procedure of the GW approximation (GWA) within the first-principles 
framework, we calculate the self energy from eigenfunctions and eigenvalues generated by the 
local-density approximation (LDA). We analyze several possible sources of error in the theory and 
its implementation, using a recently development all-electron method approach based on the full- 
potential linear muffin-tin orbital (LMTO) method. First we present some analysis of convergence 
in some quasiparticle energies with respect to the number of bands, and also their dependence on 
different basis sets within the LMTO method. We next present a new analysis of core contributions. 
Then we apply the GWA to a variety of materials systems, to test its range of validity. For simple 
sp semiconductors, GWA always underestimates bandgaps. Better agreement with experiment is 
obtained when the renormalization (Z) factor is not included, and we propose a justification for 
it. We close with some analysis of difficulties in the usual GWA procedure. 



I. INTRODUCTION 



Even though the GW approximation (GWA) of Hedini is as old as the local-density 
approximation (LDA), it is still in its early stages because of serious difficulties in its im- 
plementation. In the usual ab initio procedure, G and W are constructed from the LDA 
potential, which generate the self-energy = iG x W. Additionally the quasiparticle en- 
ergies (QPE) are usually approximated as a perturbation correction to the LDA, from the 
matrix elements of the diagonal parts of S — V^J-'^; see Eqns. ^]and IT^ below. In princi- 
ple, it is well-defined as a procedure. However, there is a controversy what the numerical 
result of this procedure is for semiconductors. Nearly always the GWA is implemented 
in conjunction with the pseudopotential (PP) approximation, which we will call PPGW. 
It was widely thought that PPGH^ predicts bandgaps in semiconductors to rather high 
accuracy. However, recent all-electron miscalculations to survey bandgaps in semiconduc- 
tors using the full-potential LMTO (FP-LMTO) by Kotani and van Schilfgaarde^ result in 
bandgaps that are generally smaller than experimental values. The result is confirmed by 
other calculations using two independently-developed full-potential linear augmented plane 
wave (FP-LAPW) codes: one by Usuda, Hamada, Kotani and van Schilfgaarde^ and another 
by Friedrich, Schindlmayr, Bliigel and Kotani^. These methods all use essentially the same 
GW codes originally developed in conjunction with the LMTO method^; they differ only 
in the input eigenfunctions. Calculations from other, independently developed all-electron 
GW methods^^ are consistent with this conclusion^. 

Tiago, Ismail-Beigi, and Louie^ used the PPGW scheme that included Si 2s and 2p cores 
in the valence to analyze the dependence of some semiconductor bandgaps on the number of 
unoccupied states N' used to construct the self-energy. They suggested that the discrepancy 
between all-electron GW and PPGVT gaps could be attributed to incomplete convergence 
in the all-electron calculations. To address this point, the convergence in A^' is taken up 
in Sec. IIIII We begin with an outline our all-electron GW method (Sec. |H|); it includes a 
comparison of the energy bands in Si to those of an APW calculation taken from Friedrich 
et al.— , establish the method's ability to reproduce near-exact LDA eigenvalues. In Sec. IIVI 
we show how selected QPEs change with increasingly larger LMTO basis sets for a variety 
of semiconductors. The results are weakly dependent on basis even for relatively small basis 
sets. We present some rationale for why this should be so, and note the implications for 
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both precision and efficiency in implementations of the GWA for basis sets in general. 

Because our results are well converged for either kind of test, we still think that PPGW 
is problematic, in contradistinction to the conclusions in a recent paper by Delaney, Garcia- 
Gonzalez, Rubio, Rinke and Godby-W, who showed that all-electron GW and PPGH^ give 
essentially the same result for the Be free atom. However, the Be atom is a special case 
in part because the all-electron and pseudo radial functions should closely correspond to 
each other (the 2p radial function has no nodes, and the only core that is orthogonalized or 
pseudized is the deep Is core (e^^^ ^ — 105 eV); moreover the PP is constructed with the 
atom itself as reference. Pseudopotentials are constructed to solve LDA reliably, but not to 
solve the GWA. There are now many detailed checks comparing PP-LDA results against 
the corresponding all-electron values; but there are few similar comparisons for GW. The 
discrepancies between all-electron and PPGl^ appear to be much smaller when PPGW 
includes the highest lying core states in the valence. 

Sec. E] analyzes different core contributions to the QPE in several semiconductors. This 
provides some insight as to what approximations may be made concerning the core; we also 
briefly consider some aspects of PPGW in this context. 

In Sec. IVII we show some new results for a variety of materials, as well as repeating some 
previously reported calculations^JiiSiHiii with rather tight tolerances. We confirm that the 
usual GWA procedure generally underestimates bandgaps. We also show that a partial self- 
consistency can be accomplished by calculating QPEs without the renormalization factor Z 
(i.e. Z=l). Semiconductor bandgaps are systematically improved using Z=l, though they 
continue to be underestimated. An important reason for this is that the LDA overestimates 
the screening of W, resulting in an underestimate of S = iGW and bandgaps. We show 
that the adequacy of the GWA varies from system to system: only when the starting LDA 
is reasonably good does the GWA reasonably predict QP energies. Thus, some kind of 
self-consistency is necessary to obtain reasonable results for wide range of materials.— 



3 



II. METHODOLOGY 



All-electron LDA in FP-LMTO 

Before turning to the analysis, we briefly describe the LDA method we use as input 
for the GW calculations. (Readers interested in conclusions of this paper not related to 
basis-set issues can skip this section). An early version of this method was presented in 



Ref. [151; we describe here how additional local orbitals are included to extend the linear 
method. Local orbitals are essential to the analysis, because QPE in GWA are sensitive to 
a wider range of states than in LDA (e.g., the LDA depends only on occupied states). One 
consequence is that the linear approximation inherent in standard linear and pseudopoten- 
tial methods is less reliable for the GWA than for the LDA. The basis functions used in 
the present technique are a generalizationi^ of the standardi^ method of linear muffin-tin 
orbitals (LMTO). Conventional LMTOs consist of atom-centered envelope functions aug- 
mented around atomic sites by a linear combination of radial wave functions (f and their 
energy derivatives ip. ip = fmi^Rhf^) is the solution of radial Schrodinger equation at site 
R at some linearization energy e^. A linear method matches the {ip, ip} pair to value and 
slope of the envelope function at each augmentation sphere boundary, which means that 
the LDA Schrodinger equation can be solved more or less exactly to first order in e — em 
inside each augmentation sphere. Envelope functions in the standard LMTO method consist 
of Hankel functions. In the present basis^^ the envelope functions are smooth, nonsingular 
generalizationsii of the Hankel functions: the /=0 smooth Hankel satisfies the equation 

{V' + e)Ho{e,r^-r) = -4ngo{r',r) (1) 
go{r'';r) = (^^/7^r'''j exp (^-(r/r*")^) ^ 5(r) as r"" ^ 

and reduces to a usual Hankel function in the limit — * 0. for higher L={l,m) are 
obtained by recursion^^. The basis can be divided into three types of functions: 

(i) A muffin-tin orbital (MTO) Xrjl which consists of a smoothed Hankel centered at 
nucleus R and augmented by linear combinations of (pm and ipm for each L channel inside 
every augmentation sphere 

X«,L(r) = HLienji, r^; r - R) + ^ C^;^^,{P«,,a'(r) - PR'k'L'{r)} (2) 

R'k'L' 



PR'k'u is a one-center expansion of Hl^erji, r\^i] r — R), and Priwu is a linear combination 
of '^Rii.^Rjh'^) and i^Bii^Rjuf) that matches Priu'v at the augmentation sphere radius. Ex- 
pansion coefficients C^^yy are chosen to make XRjL{^) smooth across each augmentation 
boundary. SRji is chosen to be at or near the center of the occupied part of that particular 
/ channel. Products of such functions enters into the construction of the hamiltonian and 
output density. The present method differs in significant ways from the usual LMTO and 
LAPW methods, and is not a simple product of MTO Eq. and bears some resemblance 
to the Projector Augmented Wave prescriptioriF^. It greatly facilitates / convergence in the 
augmentation; see Ref. 

{ii) "Floating orbitals" consisting of the same kind of function as (i), but not centered at 
a nucleus. Thus, there is no augmentation sphere where the envelope function is centered. 
There is no fundamental distinction between this kind of function and the first type, except 
that the distinction is useful when analyzing convergence. Floating orbitals make little 
difference in LDA calculations; but a basis consisting of purely atom centered envelope 
functions is not quite sufficient to precisely represent the interstitial over the wide energy 
window needed for GW calculations. Without their inclusion, errors in QPE of order 0.1 eV 
cannot be avoided, as we will show. 

{iii) A kind of "local orbital" which has a structure similar to {i). The fundamental dis- 
tinction is that the "head" (site where the envelope is centered) consists of a new radial 
function 0f evaluated at energy ef either far above or far below the linearization energy si. 
For deep, core-like orbitals 0f is integrated at the core energy; for high-lying orbitals ef is 
typically taken 1-2 Ry above the Fermi level Ep. In the former case a tail is attached, with 
its smoothing radius chosen to make the kinetic energy continuous across the "head" 
augmentation sphere. It is thus atypical of conventional local orbitals, as it is nearly an 
eigenstate of the LDA hamiltonian without requiring other basis functions. 

As is well known, the reason for using augmented wave methods in general (and especially 
the LMTO method), is that the Hilbert space of eigenfunctions in the energy range of interest 
is spanned by much fewer basis functions than with other basis sets. In the present method 
two envelope functions j = 1,2 are typically used for low / channels s, p, and d, and one for 
higher / channels (/ and sometimes g). The augmentation + local orbital procedure ensures 
that the basis is reasonably complete inside each augmentation sphere within a certain 
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energy window; the envelope functions + floating orbitals ensures completeness of the basis 
in the interstitial. The distinction between standard LMTO envelope functions and the 
smoothed ones used here is important because the generalized form significantly improves 
this convergence. Core states not treated as local orbitals are handled by integrating the 
radial Schrodinger equation inside an augmentation sphere and attaching a smoothed Hankel 
tail, allowing it to spill into the interstitial. Thus the Hartree and exchange correlation 
potentials are properly included; only the matrix element coupling core and valence states 
is neglected. 




Bands LMTO 
36-45 0.27 



26-35 0.14 



16-25 0.03 



1-14 0.01 



LAPW 
1.00 
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0.29 



0.02 
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FIG. 1: LDA energy bands in Si, computed by different methods. APW bands from Ref. 
denoted by "+" and can be regarded as near exact. Dots denote bands calculated by the same 
authors using the LAPW method, without local orbitals. Solid lines denote bands computed by the 
present generalized LMTO method, including local and floating orbitals as described in the text 
(180 energy bands were included in the basis). On the right is a table of the RMS deviation relative 
to the APW bands for several energy windows, computed along the L— P line. First column denotes 
the range of bands used to compute the RMS deviation; the horizontal dashed lines denote the 
approximate energy window for each range. Second and third columns denote the how the present 
method and the LAPW method respectively deviate from the APW bands, in eV, as described in 
the text. 



When used in conjunction with GW calculations we typically add local orbitals for states 
not spanned by {0, 0}, and whose center of gravity falls within ~ ±25 eV of the Fermi level 



6 



Ep. Both the low-lying and high- lying states can be important, and we shall return to it 
later. Fig. ^ shows the effects of linearization in Si, where an APW calculation of the LDA 
energy bands is availabloi^. Friedrich et al.— compared LDA bands generated by a full APW 
calculation to those generated by LAPW. They are reproduced in Fig. Q together with the 
bands generated by the LMTO+local+floating orbitals described above. The LAPW and 
APW are nearly indistinguishable on the scale of figure for energies up to ~25 eV. For 
energies above 25 eV, the LAPW begins to deviate from the other two, showing the effects 
of linearization. The APW and generalized LMTO bands are essentially indistinguishable 
on the scale of figure for energies up to ~40 eV. Above that, slight differences begin to 
appear; the differences gradually increase for still higher energies. Fig. Q also tabulates 
the RMS deviation from the APW bands for several energy windows. The present method 
agrees with the APW bands to ~0.01 eV for levels within Ep ± IRy, and to ~0.25 eV for 
levels below Ep + 4Ry. Friedrich et al. report similar improvements to their LAPW bands 
when local orbitals are added^. They also compare bands generated by a PP, and show 
that the errors are comparable to the conventional LAPW method. Fig. ^ establishes rather 
convincingly that the present method is nearly complete over a rather wide energy window 
in Si. When local orbitals are included, it is comparable to an LAPW method that includes 
local orbitalsi^i, and it is superior both to PP and conventional LAPW methods. 



All-electron GW with mixed basis for W 



We briefly describe our all-electron implementation of the GW approximation. A more 
detailed account will be given elsewhere^S. The self-energy E is 

S(r, r', u) = — [ duj'G (r, r', u - u') e'^^'W (r, r', u') . (3) 

In this paper G will be taken to be the one-body non-interacting Green function as computed 
by the LDA, and the screened Coulomb interaction W is computed in the random-phase 
approximation (RPA) from G. Both G and W are obtained from the LDA eigenvalues £kn 
and eigenfunctions \l/kri- For a periodic hamiltonian, we can restrict r and r' to a unit cell 
and write G as 

G.(r,r^.)^fHt.(£M. (4) 
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The infinitesimal —i5 is to be used for occupied states, and +i6 for unoccupied states. W 
is written as 

W = e-^v = {l-vliy^v (5) 

where 11 = —iG x G is the bare polarization function shown below, v = e^/|r — r'| is the 
bare coulomb interaction, and e is the dielectric function. For simplicity, the spin degree of 
freedom is omitted. 

Neglecting the off-diagonal part of S, we can evaluate QPE E'k^ from 

^kn = £kn + Zk„[(^k„|S(r,r',£kn)|^kn) - (^kn | K'^''^ W | ^kn)] • (6) 

is the quasi-particle (QP) renormalization factor 

(7) 

and accounts for the fact that S is evaluated at the LDA energy rather than at the QPE. 
Eq. (jni) is the customary way QPEs are evaluated in GW calculations. In Sec. EH we present 
an argument that using Z=l (or neglecting the Z factor) is a better choice than Eq. ((7j), 
and how the QPE is affected in actual calculations. However, the results presented here use 
the Z factor except where noted. 

In FP-LMTO, eigenfunctions of the valence states are expanded in linear combinations 
of Bloch summed MTOs, Eq. 

^kn,(r) = E«W'.^W, (8) 

RjL 

Inside augmentation sphere R, the Hilbert space of the valence eigenfunction \l/kn(i") consists 
of the pair (or triplet) of orbitals {fRi, (pRi or (pm, ip^i, ip^^i) at that site^, and can be 
represented in a compact notation {(fRu}- is a compound index for both L and one of 
the (v^ra, v^|j;) triplet. The interstitial is comprised of linear combinations of envelope 
functions consisting of smooth Hankel functions, which can be expanded in terms of plane 
wavesii. Therefore the \E'kn(i') can be written as a sum of augmentation and interstitial 
parts 

*k„(r) = E4>LW + E/?G ^G(r), (9) 

Ru G 



^kn 



d 

(^kn|^S(r,r',ekn)|^k„) 
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where the interstitial plane wave (IPW) is defined as 

P^(r) = if r G any MT 

= exp(z(k + G)r) otherwise, (10) 

and the (p\^ are Bloch sums of ipRu 

^L(r) = E^««(r-R-T)exp(zk-T). (11) 

T 

T and G are lattice translation vectors in real and reciprocal space, respectively. 

Throughout this paper, we will designate eigenfunctions constructed from MTOs as 
"VAL". Below them are the core eigenfunctions which we designate as "CORE". There 
are two fundamental distinctions between VAL and CORE: eigenfunctions: the latter are 
constructed independently by integration of the spherical part of the LDA potential, and 
they not included in the secular matrix. Second, the cores are confined to MT spherea^^. 
CORE eigenfunctions are also expanded using Eq. (jH)) in a trivial manner; /3q" = and only 
one of ^ 0. The discussion below applies to all eigenfunctions, VAL and CORE. 

Through Eq. Q, products \l/kin x ^kan' can be expanded by P^^'^^^{r) in the interstitial 
region because 

Ph\i^) X ^G.(r) = ^Gl+G.(r)- Within sphere R, wave function products can 
be expanded by B^l^^'^{r), which is the Bloch sum of the product basis {_B/jm(r)}, which in 
turn is constructed from the set of products (r) x (r)} adapting the procedure by 
Aryasetiawan^i. Eq. is equally valid in a LMTO or LAPW framework, and eigenfunctions 
from both types of methods have been used in this GW scheme^i^. We restrict ourselves to 
LMTO-derived basis functions here. 

We define the mixed basis {Mf{r)} = {P^{r) , B^j^{r)} , where the index I = {G,Rm} 
classifies the members of the basis. By construction, Mf is a good basis set for the expansion 
of products of \&kn- Complete information to calculate S and En{\i.) are matrix elements 
of the products (\l/q„|\E'q_kn'Mj ), the LDA eigenvalues e^n, the Coulomb matrix f/j(k) = 
{Mf\v\Mj-), and the overlap matrix {Mf\M^). (The overlap matrix of IPW is necessary 
because (-PgI-^g') 7^ for G 7^ G'.) The Coulomb interaction is expanded as 

^(r,r') = Y: \Mf)vu{k){M% (12) 

k,/,J 

where we define 

\Mf)^j:\Mf^)iO%], (13) 
I' 
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0),i= {M},\M}). (14) 

W and the polarization function 11 shown below are expanded in the same manner as 
Eq. jni). 

The exchange part of S is written in the mixed basis as 

BZ occ 

(^q„|S.|^q„) = EE(*q-l*q-kn'M/)^/j(k)(M^^q-kn'|^qn). (15) 
k n' 

The screened Coulomb interaction Wij{c{,uj) is calculated through Eq. (jSl), where the 
polarization function 11 is written 

BZ occ unocc 

IilM,Uj) =Y.llll W^kn|*q+kn')(*q+kn'|*knMj') 
k " n' 

X f ' ■ (I'i) 

\UJ — ^q+kn' + £^kn + T'O UJ + ^q+kn' " £^kn ~ ^0 J 

Eq. assumes time-reversal symmetry. We use the tetrahedron method for the Brillouin 
zone (BZ) summation in Eq. following Ref. |3 We first calculate the contribution to 11 
proportional to the imaginary part of the second line in Eq. ()16|) . and determine the rest of 11 
by Hilbert transformation (Kramers-Kronig relation). Such approach significantly reduces 
the computational time required to calculate 11. 
The correlation part of S is 

BZ All 

(*q„|S,(cu)|%„) = 5]5^5^(^q„|*q_wM]')(M)*q_kn'|*qn) 

k n' IJ 



oo 



x/ '-^WfAKu')^- ^ --. (17) 

Here —i5 is for occupied states; +i5 is for unoccupied states. = W — v. 
GW calculations usually approximate Eq. (jHl) by 

^kn = £kn + Zkn[(^kn|S^^^(r,r',£k„)|^kn) - (^kn | K^c^^ ( N^"^^] , | ^k„)] (18) 

where S^"^^ and V^J^^dny-^^]) are calculated only from eigenfunctions belonging to VAL. In 
the present method we calculate the Ei^n including the core contributions from 

Ekn = £kn + Zkn X [(^kn | Sx(r , r') + Sc (r , F , ^kn) | ^kn) 

-(*kn|K'c°^(K"1,r)|vI/k„,)]. (19) 
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Note V^J^^{[n^°^^^],r) is for entire density 77,*°**^^ As n' in Eq. ()15p can be divided between 
CORE and VAL, we have = E^ore ^ ^val ^j^jg p^p^^^ neglect the CORE 
contribution to Sc. In Sec. |3 we examine in some detail the contributions by shallow 
cores to correlation by including them in VAL using local orbitals, and will see that their 
contribution is small except for very shallow cores. 

In short, no important approximation is made other than the GW approximation itself; 
and it is to the best of our knowledge the only implementation of GWA that makes no 
significant approximations. Results depend slightly on what kind basis set used to generate 
G and W, as we will show, and also on the tolerances in parameters used in the GW^-specific 
part of the calculation. 

LMTO-based calculations presented here employ a widely varying set of basis functions, 
ranging from ~20— 90 orbitals per atom, as described in more detail below. They typically 
consist of a basis of spdfg+spd orbitals centered on each atom, some floating orbitals and 
sometimes local orbitals. In the Si calculations local p orbitals of either of 2p character or 
of 4p character were used, as we describe below. For the GW part of the calculation, the 
Si results shown below use parameters representative of the various system studied: LMTO 
basis functions are re-expanded in plane waves to a cutoff of 3.3 a.u. in the interstitial 
region, i.e. |k + G| < 3.3 Bohr~^ in the second term of Eq. Q. The IPW part of the mixed 
basis used to expand v, 11, and W used a cutoff |k + G| < 3.0 Bohr~^; the product basis 
part consisted of 90-110 Bloch functions/atom (we use different product basis for 
and Tj^^^). Augmentation sphere radii were chosen so that spheres approximately touched 
but did not overlap, and the product basis functions entering into the mixed basis M in 
Eqns. ()16p and (fT7|) were expanded to 1=5. In the calculation of Eq. (fT7|) the poles of G 
were Gaussian broadened by (7=0.003 or 0.01 Ry. These parameters correspond to rather 
conservative tolerances: tests at tighter tolerances in these parameters change the QPE by 
~0.01 eV. (Systematic checks were performed for each material studied.) The tetrahedron 
method was used for 11 with a 6 x 6 x 6 A;-mesh (doubling the number of points in the energy 
denominator) except where noted. The same mesh is used to calculate 11 and S. This k 
mesh is reasonably well converged, systematically overestimating conduction-band states by 
~0.02 eV in Si and similar semiconductors relative to the fully /c-converged result.— 
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III. CONVERGENCE IN QUASIPARTICLE ENERGIES: NUMBER OF UNOC- 
CUPIED STATES 




A PP (plasmon pole) 
O LMTO (converged W) 
□ LAPW 
• LMTO 



0.7 







20 40 60 80 100 120 



Number of unoccupied bands 



FIG. 2: r25u— >Xic gap in Si in GWA as a function of the number of unoccupied states A'^'. 
Filled(Yellow) squares: LAPW GWA taken from Ku and Eguilua^^. The authors presented 
data only for A^' < 31. Also, their data was given for the minimum gap, so we shifted their results 
by +0.14 eV to estimate the r25D— >Xic gap. Some checks show that the shift should be ~0.14 eV, 
approximately independent of A^'. Filled(Green) pentagons: LMTO results varying A^' both G 
and W. LMTO results were shifted by -0.02 eV to correct for incomplete k convergence^^. The 
Si 2p was included in the valence using a local orbital (the LAPW calculation of Ku and Eguiluz 
did not). The total dimension of the LMTO basis is 180. Results from the same basis are shown 
as Filled(green) pentagons in Fig. |1J and also Filled(green) pentagons No. (13) in Fig. |SJ Open 
triangles: PPGM/^ from Ref. 0, which included the Si 2p levels as part of the valence, and in which 
W was computed using the plasmon-pole approximation. Open circles: LMTO results varying 
A^' in G but not W. 

In Fig. 121 we show the r25«-^Xic gap for Si computed by Eq. (fT^ as a function of the 
number of unoccupied states A^' in Eqns. ()16p and (fT7j) : i.e. summation n' over unoccupied 
states is restricted to n' < N'. Fig. |21 depicts our main with pentagons (Si 2p treated as 
VAL). It tracks well the all-electron results of Ku and Eguilua^^, which used an LAPW 
method, except their data is ~0.05 eV less than ours. However, their results are limited 
to A^' < 31, which is not sufficient to analyze convergence for large A^'. If one assumes the 
LAPW converges with A^' at the same rate as the LMTO case, their best result with A^' = 31 
should be ~0.1 eV less than the converged value. Indeed a very recent calculation of the 
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same system by Friedrich, Schindlmayr, Bliigel and Kotani^, based on LAPW with an LDA 
basis of ~300 basis functions, showed A^'-convergence similar to LMTO. 

Tiago, Ismail-Beigi, and Louie^ presented a PPGVT calculation of some QPE in Si, Ge, 
and GaAs where they included the higher-lying core states into the valence so as to assess 
the effect of the core. They monitored the rate of convergence in QPE with A^'; their data 
for Si are shown as open triangles in Fig. |21 There are some similarities, but also two 
discrepancies: 

(i) For A^' ^ 30, the behavior is rather different. 

(ii) In the asymptotic region A^' ^ 30, the PPGiy and LMTO results converge at some- 
what different rates. 

In order to examine point (i), we tried LMTO calculations where W is fixed (i.e. A^' 
is truncated only in Eq. (jH)). This calculation (open circles in Fig. ^ tracks well the 
PP GW result for A^' ^ 30. This looks reasonable because the PP GW is combined with 
the plasmon-pole approximation, which satisfies the sum rule for Ime^^ for any A^'; thus W 
converges rather quickly with respect to A^'. However, the two LMTO calculations show little 
difference in the asymptotic behavior, which means that it is controlled by A^' in Eq. (jl7|) . 
as was already discussed by Tiago et al. 

It can be seen that the A^' dependence of either LMTO calculation is slightly different than 
the PP GW result for both intermediate and large A^': the change in the P25v—*^ic gap from 
A^' = 35 to A^' = 60 for PPGH^ is roughly twice the change obtained by the LMTO method. 
As we noted in SeclHl LMTO-LDA eigenvalues are very close to the full APW results in this 
energy range (see Fig. This indicates that the eigenf unctions are also precise. Moreover, 
Friedrich et al.— compare LDA-APW eigenvalues to an LAPW+local orbitals method; the 
three sets of eigenvalues (APW and LMTO+local orbitals, LAPW+local orbitals) are very 
close to one another. By contrast, the LDA energy bands computed by either conventional 
LAPW or PP methods correspond to APW eigenvalues far less well; see Ref. [J. Finally, the 
dependence on A^' as computed by the LAPW+local orbitals method is essentially similar to 
the present LMTO results. When these observations are considered as a whole, they suggest 
that what discrepancy does exist between LMTO+local orbitals (or LAPW+local orbitals) 
and PPGVT may be an artifact of the pseudopotential construction in the PPGW method. 
We cannot rule out possible limitations to the present method, however. Differences with 
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PPGW are small in absolute terms. Even though eigenvalues generated by LMTO and 
LAPW+local orbitals are very close to APW eigenvalues, eigenf unctions may be less well 
described. And even though the LMTO and LAPW hamiltonians are very different, the 
QPEs are generated by a common GW code. If there is some limitation in the numerical 
procedure, it would be common to both LMTO and LAPW calculations. 
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FIG. 3: T25v Xic gap in Si, as a function l/n|, where n| is the number of /c-points in the 
full Brillouin zone. The dependence on l/n| shown for T25V — *■ ^ic is essentially same for all of 
the unoccupied QPEs we examined. Gaps for nfc=6, njt=5, and rifc=4, exceed the ni.=8 case by 
0.02 eV, 0.04 eV, and 0.10 eV, respectively. We can also estimate what the — > 00 gap would be 
by extrapolating the approximately linear dependence on to zero (dashed line). The Uk = 8 
case apparently overestimates the converged result by 0.01 eV. 

It is also possible that the calculation by Tiago et al. suffers from incomplete k- 
convergence. Their PPGW used a 4 x 4 x 4 k-mesh. /c-convergence is mainly limited 
by divergent behavior for |k| ^ in Eq. ^1 To treat this divergence, we use the offset-P 
method, which was originally developed by ourselves^ and is now used by other group9^»2i. 
It is essentially equivalent to techniques that treat the divergent part analytically, as is typ- 
ically done by PPGW practitioners. Fig. El shows the dependence of P25V ~^ '^ic gap on n^, 
but in general all of conduction bands shift nearly rigidly with changes in Uk {uk = number 
of linear divisions of the /c-mesh in the Brillouin zone) . Bandgaps are approximately linear in 
the reciprocal of the total number of points, 1/nl. The figure shows that a 4 x 4 x 4 fc-mesh 
overestimates the /c-converged gap by ~0.1 eV— . This may explain most of the remaining 
discrepancy between the PPGVT calculations of Tiago et al. and the present results. 

In Sec. 13 we analyze the dependence of QPEs on the core treatment. Proper treatment 
of the core is somewhat subtle^^, and we use the local orbitals for the analysis. Because they 
are already nearly exact solutions of the LDA for the states they constructed to represent, 
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they minimally hybridize with other basis functions; consequently any higher lying CORE 
state can be readily be converted into a valence state with minimal perturbation of the LDA 
basis. Use of local orbitals enables us to investigate how different kinds of core contributions 
affect QPEs in a well-controlled and systematic way. We show that differing treatments of 
the Si 2p core only slightly affect QPEs; similar results are found for other deep cores. A 
significantly larger dependence is apparently found using the FFCVT method. 



IV. CONVERGENCE IN QUASIPARTICLE ENERGIES: BASIS DEPENDENCE 

^LDA - (eV) N' for large basis 

20 60 100 140 180 10 40 80120 160 




Number of unoccupied states N' -E'ldA ~ -^V (^^) 



FIG. 4: Left panel shows T25V gap in Si, as a function of number of unoccupied states 

N' for a smaller basis (squares) and a larger basis (pentagons). The latter are redrawn from 
Fig. 121 Top horizontal scale shows an approximate relation between energy and N' in the large 
basis (interpolated from levels at F). Right panel contains the same data but reverses the top and 
bottom horizontal scales. (Had N' in the upper horizonatal axis been drawn for the small basis 
(squares), the scale would be a little different: the last data point corresponds to A^' = 82 instead of 
120.) Inset compares convergence in Xic as a function of 1/N' to a PPGW calculation that includes 
2p states in the valence^ (triangles) and a PPGW calculation that does noli^^ (diamonds). LMTO 
data were shifted by -0.02 eV to correct for incomplete k convergence^^. The small-basis data 
(squares) in the right panel were further shifted by -0.01 eV— to clarify how large- and small-basis 
data diverge as the energy increases. 

Here we study the convergence in QPEs as the LMTO basis set changes, retaining all the 
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eigenfunctions for a given basis in the calculation (A^' encompasses all unoccupied states). 
A given LMTO basis defines a finite Hilbert space of eigenfunctions; the GWA is a well- 
defined procedure in that space, and we can study how the QPEs change as the Hilbert space 
is refined. This procedure corresponds more closely to analyses of basis set convergence 
common in other kinds of calculations (e.g. LDA and Hartree-Fock). We can also anticipate 
that it will be smoother than the A^' truncation of Sec lIIH indeed this will turn out to be the 
case (see especially Fig. El): the band gaps are insensitive to the choice of basis once a certain 
level of completeness is reached. It its also obviously true that the Hilbert space depends 
on the choice of basis constructing it. Therefore, the results presented here are specific to 
the LMTO basis described in Sec. IHl and in particular what kinds of orbitals are included, 
e.g. orbitals oi f ot g character, or local orbitals that correct the linearization inherent in 
most of the standard methods (LMTO, LAPW, and the construction of a norm-conserving 
PP). By adding different kinds of orbitals we can identify how different parts of the Hilbert 
space, (most notably corrections to linearization common to most methods) affect QPEs. 
Since the LMTO basis is tailored to the crystal potential, the LDA eigenfunctions converge 
more rapidly with the basis dimension than do plane-wave based basis sets^i. Consequently 
we might expect a more rapid convergence in the GWA QPEs. On the other hand, by 
transformation to, e.g., Wannier functions, it should be possible to design a generic scheme 
that exhibits similarly rapid convergence. 

Initially, we compare in Fig. |3 the dependence of the r25t,^Xic gap in Si on A^' for 
two basis sets: one relatively small and another relatively large. The right panel of Fig. |3] 
compares the same data, but the horizontal scale corresponds the energy of state A^' at 
r. Also in this panel, the small basis-data at small energy was artificially shifted down 
by 0.01 eV to make it easier to see how the two cases depend on energy.—. The difference, 
initially 0.01 eV, increases by an additional 0.01 eV as the energy E approaches ~ Ef +50 eV. 
For higher energies, the discrepancy between the two increases more rapidly. This is because 
the ability of the small LMTO basis to describe eigenvalues above Ef + 40 eV begins to 
degrade. 

However, we can see that the gaps at the respective maximal A^' (i.e. all unoccupied 
states included) are in good agreement. Including all unoccupied states in a limited basis 
is a another kind of hilbert space truncation, but it is also well defined: use use the Hilbert 
space of eigenfunctions basis consisting of LMTO eigenfunctions and their products (Sec. lIT|) . 
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and apply the GWA within that Hilbert space. The LMTO basis can efficiently choose the 
important part of the Hilbert space tailored to the crystal potential. Thus good agreement 
need not be some fortutitous artifact of this particular pair of the LMTO basis sets, even 
though the maximal A^' is small in light of a traditional N' cutoff analysis22,. Indeed the 
A^' cutoff of Sec mil may choose the Hilbert space less well, especially since that kind of 
truncation is not smooth. Below we present detailed analysis of the dependence on the basis 
set to justify the good agreement in Fig. the band gaps are insensitive to the choice of 
basis once a certain level of completeness is reached. 

3.3 r - 
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FIG. 5: QPE in the GWA at Fisc (top), Lie (middle), and Xic (bottom) in Si relative to the 
valence band maximum, using different basis sets in the present FP-LMTO. Abscissa is the total 
number of basis functions A^. Yellow diamonds are a minimal basis (see text). All results depicted 
by fill circles contain no local orbitals. Those depicted by filled (empty) squares or pentagons 
include the Si 2p (Si 3p) as local orbitals. See text for further description. 

Fig. El shows results of a systematic study of the convergence in the first unoccupied QPE 
at F, X, and L in Si with progressively larger basis sets. LDA eigenvalues are not shown 
because they are the same within ~ 0.01 eV for all cases (0.60 eV for X, 1.42 eV for L, 
2.52 eV for F). These data comprise very diverse basis sets, particularly for the LMTO 
method which traditionally uses a minimal basis. Some details concerning these sets help 
explain in what manner convergence is reached: 



]1 Filled (yellow) diamonds (1) includes spdfsp atom-centered functions, and is the 
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only basis without floating orbitals. There are no local orbitals; Si ls,2s,2p are CORE. 

§2 Filled(blue) circles (2) and (3) add floating orbitals of sp and of spd character, 
respectively. Their effect is to cause QPEs to decrease slightly relative to (1). Adding still 
more floating orbitals (even large numbers of them) shift QPEs by ~0.01 eV. 

§3 Other Filled(blue) circles (5, 6, 10, 12) include still more envelope functions 
comprised of a mixture of atom-centered functions and floating orbitals, but adding no 
local orbitals. 

§4 Filled(green) Square and Open Square (4, 7) correspond to (3, 6), respectively, 
but adding a local orbital (Green: Si 2p, Open: Si 4p). When the 2p is included as VAL, 
CORE consists of Sils,2s only. A local orbital of either Si2p or Si4p shifts QPEs — in 
roughly equal but opposite directions. 

§5 Filled(green) Pentagon and Open Pentagon, (8,9,11,13) include an additional 
Si4(i local orbital. (11) corresponds to (10) + (Si2p or Si4p)+ Si4(i. (13) corresponds to 
(12)+... as well. (8) is (6)+Si4p+Si4d. The effect of Si4d is small. 

These points show in a compelling way that once the basis reaches a certain level of 
completeness, the change of QPE with further enlargement is very small. Set (1), which 
consists only of atom-centered functions, is somewhat incomplete except inside the aug- 
mentation spheres where the eigenfunctions are constructed out of linear combinations of 
{0,0}. Considering the open structure of zincblende, such a basis may be expected to be 
less complete in the interstitial. Comparison basis sets without local orbitals (circles) with 
set (1) shows that this particular purely atom-centered basis is slightly deflcient for reliable 
calculation of QPE, since the addition of floating orbitals induces a (/c-dependent) reduction 
in the conduction band of ~ 0.02 — 0.10 eV. It is an open question whether a still more 
sophisticated atom-centered basis^^ would be adequate to describe the interstitial. 

Once the interstitial is reasonably complete (c/. sets (3) and higher), there is an almost 
negligible dependence on basis provided no orbitals are included that extend the linear 
method or alter how the core is treated. Basis sets marked by a common symbol (squares, 
circles, pentagons) share essentially the same Hilbert space in the augmentation regions; only 
the basis set corresponding to the interstitial region changes. The variation is ±0.01 eV for 
a wide range of basis sets^S. 
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Fig. El also gives us some insight to the hmitations of the hnear method. Basis sets 
(3) and (4), which differ only in how the Si p channel is treated inside the augmentation 
region, affect QPEs more strongly than radically enlarging the Hilbert space of the envelope 
functions — compare (3)^(4) and (3)^(12). Envelope functions affect only the interstitial; 
they negligibly affect the Hilbert space of the augmentation region. For the latter it is 
largely irrelevant how many envelope functions are used — and consequently the size of N' 
entering into Eqns. ()16|) and (fT7|) . What is relevant is the completeness of {</>, <p}, and results 
are independent of basis dimension provided that the entire {0, 0} Hilbert space is included. 
Said another way, the LMTO method is by design reasonably complete over a certain energy 
window in the augmentation spheres, more or less independent of the envelope functions. 
A similar story may be told for the interstitial: sets (3, 6, 10, 12) differ in the number of 
envelope functions by as much as a factor of three, but QPEs are unchanged within ±0.01 eV. 
QPEs shift when {0, 0}^{0, 0, 0^^} or {0, 0}— >{0, 0, 0^^}, essentially independently of the 
number of envelope functions (3^4, 6^7, 10^11, 12^13). 

Table HI shows three other materials (CdO, Ge, and GaAs). We can see (1) rapid conver- 
gence in QPEs as the basis is enlarged for a fixed set of augmentation functions; and (2) 
extensions to a linear augmentation affect QPEs in an manner approximately independent 
of the total dimension of the Hamiltonian. (In GaAs, both 3d and Ad must be included 
as VAL. If not, significant errors resullj^). We have tested a number of other materials as 
well, and these trends appear to be rather general. As might be expected, the number of 
basis functions needed to make the Hilbert space reasonably complete depends somewhat 
on the elements involved. The heavier atoms have larger radii and consequently slower 
/-convergence in the number of envelope functions needed; also d orbitals often play an im- 
portant role. More orbitals are required to make the basis complete when heavier atoms are 
involved. 

As noted, the linear {0, 0} Hilbert space is already reasonably complete in the case of 
Si. But this is not true in general: oxides and nitrides form a materials class where the 
effect is significantly larger. CdO is one such example (CdO forms in the NaCl structure; 
the valence band maximum falls at L and the conduction band minimum falls at F.) As 
happens for Si, there is a weak dependence on basis when the number of envelope functions 
is changed and the Hilbert space of the augmentation is held constant. But the QPEs change 
by ~ 0.15 ± 0.05 eV when the O 3p and Cd 5d states are added, as Table H] shows. (In this 
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TABLE I: QPEs of the first unoccupied state at F, L, and X, for different basis sets, in eV (relative 
to valence maximum). Columns Ua, nf, and ni denote the number of atom-centered functions, 
the number of floating orbitals, and the number of local orbitals respectively. The hamiltonian 
dimension is the sum of these numbers. Experimental data are adjusted for spin-orbit coupling by 
adding 1/3 of the splitting in the Fis valence bands. The first four CdO basis sets are identical to 
the last four except for the addition of local orbitals in the O 3p and Cd 5d channels. A 6 x 6 x 6 
A:-mesh was used in these calculations. 
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particular case it is the O 3p contribution that is dominant; however, cases arise when the 
contributions from high- lying d oi f orbitals can be of order 1 — 2 eV. NiO is such a casein.) 

The inset of Fig. E] shows some PPGW results for reference. Based on the observation 
that cutoff in the Hilbert space should be important, PPGl^ data by Tiago et al. should 
be extrapolated to 1/A^' — because they used a very large LDA basis. 



V. CORE CONTRIBUTIONS TO Sc 



In a series of papers, Shirley and co-wor 



cers analyzed the effects of the core on QPE 



in atoms (Shirley, Mitas and Martin, Ref. |3J, Shirley and Martin, Ref. |35|) semiconductors 



(Shirley, Zhu and S. G. Louie^SiS) within the pseudopotential framework. Approximate 
core contributions to both Eq. (fT^ and p7|) were evaluated. They also compared pseu- 
dopotentials constructed from both LDA exchange and from Hartree-Fock (HF) exchange 
for atoms and molecules^, and incorporated pseudopotentials of both types in studying core 
effectsS^iSi. They found sizable shifts in QPE in Si, and rather dramatic and /c-dependent 
shifts in Ge and GaAs. These analyses highlight the importance of core effects. However, 
the decomposition of the various core 

is rather closely tied to the pseudopotential construction that was a part of their implemen- 
tation. This makes it a little difficult to disentangle the various contributions. 

Here we examine contributions from the shallowest cores to Sc within the framework of 
our GW. As we noted, all the eigenf unctions are divided into two groups CORE and VAL 
as explained in Sec. IHI Using local orbitals we can represent the shallowest cores in VAL. 
To distinguish true core effects from artifacts of implementation^^, we include these cores 
in the valence with local orbitals and treat them in a special way, as described below. We 
denote such eigenf unctions as core, and the rest as val. Thus we distinguish three kinds of 
orbitals: CORE, core, and val: 

All eigenfunctions 

(20) 

CORE VAL 

core val 



In Si, for example, we use: CORE={ls, 2s}, core={2p} and val={3s,3p,3d, ...}. Because 
the core states are well separated from higher lying states, G can be partitioned into G = 
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TABLE II: QPE of the first unoccupied state at T, L, and X relative to top of valence, for core 
treatments (i) — (iv) as described in the text, in eV. States and corresponding eigenvalues e^^"^ 
treated as core are: 2p in Si, 3d in Ga and Ge, and 2s in Mg. Si data corresponds to basis set (13) 
in Fig. El GaAs data corresponds to the 68+18+6-orbital basis in Tabled Ge data corresponds 
to the 68+18+10-orbital basis in that table. Here G means [G""'''' + G""-^), W means W[U] (see 
text). A 6 X 6 X 6 fc-mesh was used in these calculations. For results with better /c-convergence 
and larger basis sets, see Table HTll 
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qCOKE _|_ Qcore _|_ Qval ^ always Calculated from the entire G, while Sc is calculated 

from core and val only (we do not consider any case where some portion of the self-energy is 
supplied by the LDA): Sc = i^G"""'^ + G™')^^ where = W[U] - v, and U is calculated 
from G'^°^'^ + G^^'K Thus core states contribute to Ec directly through G in iGW^, and also 
through W^. We resolve these contributions; that is, we calculate Sc in one of four ways: 

(i) neglect the core contribution to Sc entirely: i.e. Ec = iG'"^'(Vr[n''"'] — v), where 
iy[n'"*'] means that 11 is calculated from G""' only. We denote this as "exchange-only 
core" . 

(ii) neglect the core contribution to screening : Sc = i[G'^°^'^ + G™')(iy[n™'] — v). 

[iii) neglect the core contribution to G: E^. = zG™'(W^[n] — v). 

(iv) Sc = i{G'^°^'^ + G™')(Vr[n] — v): there is no distinction between core and val states. 

Table ITTlshows that the difference between exchange-only (i) and GW (iv) approximations 
to core treatment is small in Si (~0.03 eV for Xic). As expected, the adequacy of an 
exchange-only core depends on how deep the core is. The exchange-only approximation for 
shallow cores, such as the Ga 3d and In Ad, and the highest lying p core in the column I 
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(Na, K, Rb) and column IIA alkali metals (Mg, Ca, and Sr), is rather crude. It is interesting 
that the core contributions to 11 and to G are not always additive. 

The difference between {Hi) and {iv) is in general rather small; that is, inclusion of core 
contributions of 11 alone is sufficient to bring QPE results within 0.05 eV of the full results 
in Table im except for the very shallow Ga "id channel. For moderately deep cores, exchange- 
only treatment (i) is generally adequate, as Aryasetiawan suggested. A rough rule of thumb 
seems to be: for cores whose total charge Qspiu outside the augmentation radius is less than 
0.01 electrons, exchange-only treatment of them results in errors ~0.1 eV or less for the 
lowest excited states. (This radius may be taken as approximately half the nearest-neighbor 
bond length). 

Inclusion of core contributions to 11 can significantly increase the computational cost (in 
the Si case, leaving out the 2p contribution to 11 reduces the computational cost by ~40%). 
The relative smallness of corrections to exchange-only treatment, and the observation that 
core contributions to 11 alone are adequate for all but the most shallow cores, suggests that 
a simple approximate inclusion of core contributions to 11, Eq. (fT?)|) . should be adequate 
for all but the most shallow cores such as the Ga "id. (Fleszar and Hanke proposed a 
construction for pseudopotentials when core states are not pseudized^^.) Supposing the 
core were confined to the augmentation sphere at site i?, we can eliminate all contributions 
to the matrix element (M/'^l/knl^q+kn') except from the product-basis contribution at R. 
Since also the augmented part of \1/ depends rather weakly on ipm we can neglect the ipm 
contribution to the eigenfunctions and assume that (Mj \l/kn|^q+kn') only depends n, n', 
k, or k + q through the coefficients, {oi^^*a^^^ . Moreover, the core level energy is large 
and negative, and nearly independent of k or n. Since the dominant contributions to 11 will 
come from coupling to low-lying states, we can approximate Skn — ^k'n' by a constant, e.g. 
Ep — ecorc- These approximations are all modest but can vastly simplify the computation of 

j];CORE 

The fact that the core spills out slightly from the augmentation region needs to be taken 
into accounli^^. This can readily be accomplished by integrating the core and corresponding 
valence ^pi to a larger radius, and orthogonalizing ipi to the core. Checks show that the 
adjustment to 0/ is small unless the core is very shallow, in which case the core should be 
treated clS cl valence state. 
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Comparison to PPGW 

When the highest cores are put exphcitly into the valence as Tiago et al. did, there is 
reasonable agreement between PPGW and our results for sp semiconductors. Comparison 
with the paper of Tiago et al. to Table IH] below shows that there is agreement at the 
~0.1 eV level in Si^ and similar agreement is found for GaAs and Ge, with the PPGl^ 
results systematically higher than our results by ~0— 0.1 eV. Similarly, Fleszar and Hanke^^ 
calculated QPEs in the GWA for a variety of II- VI semiconductors, including the highest 
s and p cores in the valence. Their values are also in reasonable agreement with results 
presented in Table ITTTl fbelow) . though the PPGW data are systematically higher by ~0.0- 
0.2 eV. (Part of the discrepancy can be traced to contributions from high-lying d states, 
which are included in the present calculation using local orbitals.) Even when the high- 
lying s and p core states are included explicitly in the valence, it still seems to be the case 
that PPGiy band gaps are systematically slightly larger (by ~0.1 eV) in semiconductors 
than our GW predict. 

Materials involving transitions metals are rather more complicated. In a recent PPGW 
calculation, Marini, Onida and del Sole analyzed the QP valence bands of Cv^, comparing 
in some detail the occupied d bands to photoemission experiments. The LDA places the 
position of these levels approximately 0.5 eV closer to Ep than the experiments show. The 
authors find that the d bands narrow and shift downward by approximately 0.5 eV, bringing 
the PPGW d bands into excellent agreement with photoemission experiments. They report 
that the PPGW results depend rather dramatically on the treatment of the Cu core 3s 
and 3p levels: that it is necessary to include both states explicitly in the valence to obtain 
reasonable results. They found that the correlation contribution SJ;™'^ from these states 
shifts the d bands downward ~0.5 eV. 

We conducted a similar calculation using the present all-electron GW, and find a very 
different result. In our case, the GW correction to the LDA d bands is small — between and 
0.1 eV. Moreover, QPEs are essentially independent of how the Cu 3p state is treated: the 3d 
levels change by less than 0.05 eV when the Cu 3p state is explicitly included in the valence 
(using a 3p local orbital), as compared to being treated as core at the exchange-only level. 
The Cu 3p state is rather deep, and the weak dependence on correlation contributions from it 
is consistent with the rule of thumb indicated above: Qspiu ~ 0.005 electrons; £3^"^ ~ —70 eV. 
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In the Cu case, it appears likely that the main discrepancy between PPGW and our GW 
(whether d bands shift by 0.5 eV or not) originates in the discrepancies in S^°''°. 

VI. ADEQUACY OF GWA APPLIED TO A RANGE OF MATERIALS 

In Ref. Q Delaney et al. argued that GWA based on the LDA eigenfunctions and eigen- 
values, is an adequate (or better) approximation than self-consistent GW. It is apparently 
the case that self-consistency worsens agreement with experiment for the Be atom. More- 
over, Holm and von Bartb^ found that the valence bandwidth of the homogeneous electron 
gas is considerably worsened by self-consistency; similarly a self-consistent GW calculation 
Ku and Eguiluz resulted in an overestimate of the valence bandwidth in Ge^. Thus, self- 
consistency of this type has shortcomings. On the other hand, even in simple materials such 
as sp semiconductors, GWA bandgaps based on LDA eigenvalues and eigenfunctions are 
always underestimated when properly calculated^*^. The GWA based on LDA is evaluated 
as a perturbation relative to LDA; thus the band gap can be poor if the LDA itself is poor. 
Thus, some kind of self-consistency is necessary to reduce the dependence on the starting 
point. 

In this section, we consider three points about the GWA based on LDA, Eq. 

(A) Use of the Z factor. We show that using Z=l in Eq. (jHl) is a way to include partial 
self-consistency, and it should be a better approximation than including the Z factor. 

(B) Off-diagonal S. Eq. © is a perturbation treatment that involves only the diagonal 
matrix element of S. We consider the effect of the full S in a variety of systems analyzing 
how the adequacy of GWA is dependent on the adequacy of LDA. Even GWA with Z=l 
fails for cases when the starting LDA is poor. 

(C) Band-disentanglement problem. Even when LDA eigenfunctions are reasonable, 
if eigenvalues are wrongly ordered the perturbation treatment can have important adverse 
consequences. 
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A. Z factor 



Let us consider a partial kind self-consistency where only eigenvalues are updated: both 
eigenfunctions and W are unchanged from the LDA. This is a little different from the 
customary usual eigenvalue-only self-consistency, where eigenfunctions are frozen but W is 
updated. Updating eigenvalues widens semiconductor bandgaps. This reduces the screening, 
which causes W to increase, which in turn cause gaps to increase still more. Thus we expect 
that results from such the kind of partial self-consistency we are considering here should fall 
somewhere between the usual one-shot GW and the usual eigenvalue-only self-consistency. 
Partial self-consistency, while incomplete, should result in better QPEs than the standard 
1-shot GW ^ since eigenvalues shift in the right direction. Appendix 1X1 evaluates how Eq. © 
gets modified for a model two-level system. The result is that this kind of self-consistency 
can be approximately realized by putting Zk„ = 1 in Eq. (jH)). A different justification 
for omitting the Z factor emerged from a paper of Niquet and Gonze^, who calculated 
the interacting bandgap energy (within RPA) to obtain a correction to the Kohn-Sham 
gap. They found that the difference is essentially Eq. © with Z=\. Finally, a further 
justification for using Z=\ is discussed in chapter 7 of Ref. |44. Z=\ corresponds to the 
Rayleigh-Schrodinger perturbation, Z with Eq. (??) to the Brillouin-Wigner perturbation. 
It shows the Z=\ scheme should be better for the Frolich Hamiltonian Mahan analyzed. 

The calculations in Table UTTl support the argument that using Z=\ is a better approx- 
imation than including Z: semiconductor bandgaps are in significantly better agreement 
with experiment. They continue to be smaller than experimental values, which can be qual- 
itatively understood as follows. Using Z=\ corresponds to updating G, but leaving W de- 
termined from the LDA eigenfunctions and eigenvalues. Because the gap is underestimated 
in the construction of 11 and W ^ 11 is overestimated so that W is screened too strongly; thus 
S is too small. It is interesting, however, that QPEs evaluated with Z=\ can be rather good 
at times because of a fortuitous cancellation of errors. We can refine self-consistency by up- 
dating in a manner similar to the updating of G: that is using eigenvalues from Eq. (jHl) in 
the calculation of 11. However, e computed in the RPA (e = 1 — fll), omits excitonic effects. 
Inclusion of electron-hole correlations to 11 (via e.g., ladder diagrams) increases Ime(ci;) for 
LO in the vicinity of the gap for semiconductors. There is a concomitant increase in Re eiiS) 
for cij — > 0, as is evident by the Kramers-Kronig relations; see e.g., Ref.Q. Errors resulting 
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TABLE III: Fundamental gap, in eV. (For Gd, QPE correspond to the position of the majority 
and minority / levels relative to Ep; for Cu QPE corresponds to the d level.) Low temperature 
experimental data were used when available. QPEs in the "GVF" column are calculated with 
usual GWA Eq. © and Eq. ((7j). In the "Z=l" column the Z factor is taken to be unity. In the 
column the off-diagonal parts of S are included in addition to taking Z=l. /c-meshes of 
8 X 8 X 8 A; and 6x6x6 were used for cubic and hexagonal structures, respectively (symbol w 
indicates the wurtzite structure). GW calculations leave out spin-orbit coupling and zero-point 
motion effects. The former is determined from A/3, where A is the spin splitting of the Ti^y level 
(in the zincblende structure); it is shown in the "A/3" column. Contributions to zero-point motion 
are estimated from Table 2 in Ref. |13 and are shown in the "ZP" column. The "adjusted" gap 
adds these columns to the true gap, and is the appropriate quantity to compare to GW. 
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from the neglect of excitonic contributions to e partially cancel errors resulting from LDA 
eigenvalues, as shown by Arnaud and Alouani^. Thus, W calculated from LDA eigenvalues 
is not so bad in many cases because of this cancellation. Often calculated from the LDA 
eigenvalues is better than e^o calculated from LDA eigenvalues shifted by a scissors operator 
to match the experimental band gap (See TABLE III in Ref. |46|). This cancellation means 
that GW{Z=1) can be often be rather good, since W itself is also better than what would 
be obtained from (eigenvalue) self-consistency. Table IIIII shows that the fundamental gap 
for GW{Z = 1) is quite good for mostly-covalent semiconductors such as Si or GaAs, but 
that the agreement deteriorates as the ionicity increases. 



B. Off-diagonal contributions of S 

The usual GWA in Eq. does not include the off-diagonal contribution of S — V^^^^. A 
simple way to take into account the contribution off-diagonal parts is to replace the energy- 
dependent matrix S with some static hermitian matrix V^^ as in the following, and to solve 
the eigenvalue problem, replacing V^^^ in the LDA hamiltonian with this potential. We 
take 

V^^' = J E 1^^) + Re[S(£,)]..} (^,1, (21) 

for S. Here Re signifies the hermitian part, the eigenvalues Si and the eigenfunctions ipi are 
in LDA. This V^'^ is used in our QP self-consistent GW methodii^i^J^. This V^^ retains the 
diagonal part contribution as in Eq. © (we now consider Z=l case). From the perspective 
of the QP self-consistent GW method, including the off-diagonal S corresponds to the first 
iteration, and the LDA corresponds to the 0th iteration. Table lllll shows how the funda- 
mental gap is affected by the off-diagonal parts of V^'^ for selected semiconductors. Because 
the semiconductor eigenfunctions and density are already rather good, the off-diagonal con- 
tributions are small. Contributions from the off-diagonal part of V^'^ significantly increase 
when eigenfunctions have significant d character (see SrTiOa and ScN in Table IIIIll . For 
correlated systems the effects can be rather dramatic; see Ref. 112 how the QPEs are affected 
by the off-diagonal parts of S in Ce02. 

In general, GWA errors are rather closely tied to the quality of LDA starting point. 
In the covalent sp semiconductors C, Si and Ge, GW gaps are rather good for Z=l. In 
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the series Zn(Te,Se,S,0), the deviation between the LDA and experimental gap steadily 
worsens, and so does the GW gap. For ZnO, and even more so in CuBr, the GW gap falls 
far below experiment. For these simple sp materials, errors are related to their ionicity, 
which can be seen qualitatively as follows. As ionicity increases, the dielectric response 
becomes smaller; consequently the nonlocality missing from the LDA exchange-correlation 
potential^ becomes progressively more important. Roughly speaking, a reasonable picture 
of electronic structure in sp systems resembles an interpolation between the LDA, which 
has no nonlocality in the exchange and underestimates gaps, and Hartree-Fock, which has 
nonlocality but wildly overestimates gaps because the nonlocal exchange is not screened. As 
ionicity increases the gap widens and the dielectric function decreases. As the screening is 
reduced the LDA becomes a progressively worse approximation. Thus, the LDA is not an 
adequate starting point for GW in the latter cases. 

Discrepancies between GW and experiments become drastic when electronic correlations 
are strong. The GWA band gap for the antiferromagnetic-II NiO is far from experiment, 
and moreover the conduction band minimum falls at the wrong place (between F and X— ). 
As Table ITTTl shows, the LDA puts / levels in Gd too close to Ep- GWA results are only 
moderately better: shifts in the Gd / level relative to the LDA are severely underestimated 
(see Table lln|l . 

GW based on the LDA fails even qualitatively in CoO: it predicts a metal with Ep passing 
through an itinerant band of d character. In this case, the GWA gives essentially meaningless 
results. To get reasonable results it is essential to apply the GWA with a starting point that 
already has a gap. To get band gap for CoO in the single band picture, a non-local potential 
which breaks time-reversal symmetry is required, something which is not built into the local 
potential of the LDA. Similar problems occur with ErAs: the LDA predicts a very narrow 
minority / band straddling Ep, whereas in reality the minority / manifold is exchange-split 
into several distinct levels well removed from Ep^. GWA shifts the minority / levels only 
slightly relative to the LDA: the entire / manifold remains clustered in a narrow band at 
the Fermi level, appearing once again qualitatively similar to the LDA. 

Generally speaking, the GWA (even GWA{X=1)) is reasonable only under limited 
circumstances — when the LDA itself is already reasonable. 
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C. Band disentanglement problem 
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FIG. 6: Energy bands in Ge for k = 27r/a [00k] for small k within the GWA, using Z=l. Spin 
orbit coupling was omitted. Three approximations are compared: LDA (dashed blue line), GW in 
the diagonal-S only approximation, Eq. © (black line with circles), and GW with T, computed 
according to Eq. (|21() (solid green line). In all three cases the three states of p character (r25' 
symmetry) form the valence band maximum; this was taken to be the energy zero. The LDA 
predicts the conduction band, the T2' state of s character, to be slightly negative, causing the 
energy bands to be wrongly ordered at T. For k > 0, the T25' state of pz symmetry couples to 
the T2' state, and the two repel each other. Both kinds of GW put T2'c at approximately the 
correct position, 1 eV. However, the diagonal-only GW must follow the topology of the LDA: 
the eigenvectors are unchanged from the LDA. Therefore the band starting out at T2'c sweeps 
downward, while the Pz band starting out at T25'v sweeps upward, and the two bands cross near 
k = 0.02. When the off-diagonal parts of T, are included, these two bands repel each other as they 
should. 

Even for the simple sp semiconductors, there can be a "band disentanglement problem" 
as a consequence of the diagonal-only approximation. At times the LDA orders energy levels 
wrongly: in hep Co, for example, it inverts the order of the minority Fs and F3 levels, which 
correspond to states of L3 and L2' symmetry in the fee structure. Wrong ordering of levels 
is a particularly serious difficulty for narrow-gap semiconductors such as Ge, InAs, InSb, 
and InN. Because the LDA underestimates bandgaps, the energy band structure around F 
has an inverted structure: the s-like conduction band of Fi symmetry (labeled as F2' in the 
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homopolar case) incorrectly falls below the p-like states of r25' symmetry. 

When the GWA is evaluated from Eq. ©, the energy bands retain the same connectivity 
as in the LDA, as Fig.Elshows. Consequently the conduction band has a nonsensical negative 
mass near F, and it crosses with one of the valence bands. The diagonal-only approximation 
cannot make Ge an insulator in principle, even though the levels are properly ordered at 
r. This problem is avoided if the off-diagonal parts of S are included, as Fig. IHl shows. 
The conduction-band effective mass in the latter case is computed to be m* = 0.042mo, 
in good agreement with a value of m* = 0.038mo estimated from magnetopiezorefiectance 
spectr This shows that the off-diagonal contributions of S are reasonably well described 
by Eq. 

VII. CONCLUSIONS 

To conclude, we have analyzed various possible sources of error in implementations of the 
GWA, using calculations based on an all-electron method with generalized Linear Muffin Tin 
Orbitals as a basis. We analyzed convergence in QPEs with the number unoccupied states 
A^': the rate of convergence for intermediate A^' (where the LMTO energy bands were shown 
to precisely replicate APW bands), was qualitatively similar to, but roughly twice that of a 
PP analysis by Tiago, Ismail-Beigi, and Louie^. On the other hand, it closely tracked the 
convergence calculated by an LAPW+local orbitals method, which had a very similar LDA 
band structure. More generally, our GWA that properly subtracts V^J^^ calculated from 
the full density are in reasonable agreement with each otheri^iMJi^; those that subtract 
valence density only^ are also in reasonable agreement for cases such as Si and SiC where 
the cores are sufficiently deep. Our own experience suggests that the LDA treatment of 
core levels, where QPE are computed from Eq. (jl8p . will be problematic for GW— unless 
the cores are very deep. Since a PP construction is an approximation whose justification 
is grounded in an all-electron theory, we should expect GW calculations based on an LDA 
PP should be similarly problematic. There is apparently a significant dependence on how 
cores are treated in PP implementation o^^i^^i^^i^^ , even in Si and Cu with their deep 2p and 
3p cores. 

We then presented a new analysis of convergence that is of particular importance for 
minimal-basis implementations, and argued that measuring convergence in the traditional 
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cutoff procedure — by number of unoccupied states A^' as given in Figs. |21 and E] — are not 
particularly meaningful for a minimal basis. We presented an alternative truncation of the 
full Hilbert space of eigenfunctions, namely to use the entire hilbert space of a relatively small 
basis. We showed that a suitably constructed minimal basis is sufficient to precisely describe 
the GWA QPE within 1 Ry or so of the Fermi level, and that this kind of cutoff procedure 
seems to be more efficient than the traditional A^' cutoff of a large basis. We also showed 
that traditional linearization of basis functions, either explicit in an all-electron method 
or implicit through the construction of a pseudopotential, result in errors approximately 
independent of the size of basis. Addition of local orbitals to extend the linear approximation 
results in modest shifts in sp nitride and oxide compounds, and shifts of order 1-2 eV in 
transition-metal oxides. 

We analyzed core contributions to the self-energy, and showed that an exchange-only 
treatment of the core is adequate in most cases. For all but the most shallow cores (such 
as Na 2p and Ga 3d), we showed that it is sufficient to include the core contribution to the 
polarization only; an approximate and rather painless implementation was suggested. These 
results can provide a framework for improved treatment of the core within a pseudopotential 
approximation. 

Finally, we considered the adequacy of GWA based on the LDA, for different kinds of 
materials, and also Eq. © as an approximation to the GWA. We presented logical and 
numerical justifications that using Z=l generally gives better band gaps in insulators. In 
general inclusion of the off-diagonal part of S and some kind of self-consistency is essential 
to make the GWA a universally applicable and predictive tool. Taking into account both 
theoretical and practical aspects, the quasi-particle self-consistent GW scheme we have 
proposedii*i2*i^ has the potential to be an excellent candidate for such a tool: it obviates 
some of the difficulties seen in the standard self-consistency, it no longer depends on the 
LDA, and it appears to predict QPEs in a consistently reliable way for broad classes of 
materials. 

This work was supported by ONR contract N00014-02-1-1025. S. F. was supported 
by DOE, EES Contract No. DE-AC04-94AL85000. We thank the Ira A. Fulton High 
Performance Computing Initiative for use of their computer facilities. 
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APPENDIX A: JUSTIFICATION FOR Z = 1. 



Let us consider a limited self-consistency within GWA as follows. We restrict self- 
consistency as follows: 

1. We make only the QPE self-consistent. Eigenf unctions are constrained to be the LDA 
eigenf unctions. 

2. W is assumed to be fixed. Thus only the eigenvalues entering into G are made self- 
consistent. 

Under these assumptions, we can show that QPE are rather well approximated by Eq. (jH)) 
with Z=l. To illustrate it, consider a two-states model whose LDA eigenvalues and eigen- 
functions are given by ipi,ei, and ip2,£2, and the Fermi energy falls between these states: 
62 > E-p > El. Then the LDA Green's function is 

= -MM^ + JtllM^, (Ai) 

Lu — El — id Lu — E2 + io 

After the limited self-consistency is attained, we will have eigenvalues: 



Giuj) = + (A2) 



11 

uj — El — i6 ' uj — E2 + iS' 



where Ei is given by 

Ei=Ei + Re(^i|S(Ei, [G]) - Kc°^|^i). (A3) 

There is a similar equation for E2. Note that T,{Ei, [G]) is calculated in GWA from G of 
Eq. (HSl) at El. 

As we can expect that W is dominated by diagonal terms Wi{uj) = {ipiipilW {uJ)\^pl'^pl) and 
W2{uj) = {ip2'^2\W{uj)\ip2'ip2) , we neglect other matrix elements of W{uj). Then S becomes 



Re(^i|(S(Ei,[G])|^i) = Re I {iJi\iG{Ei + u')W{u')\^i)duj' 

iWi(uj') duj' 



Re 



Re 



El + uo' - El - i5 
iWi{u')duj' 



El + Uj' — El — i6 

Re(z^i|s(ei,[GL°^]) (A4) 
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A similar equation applies for i?2- The energy shift Ei — > ei entering into the evaluation S 
is exactly compensated by the energy shift in G ^ G^^^. Or equivalently using Z=l is an 
approximate way to obtain self-consistency. Eq. ()A4|) corresponds Eq. (jH)) with Z=l. 
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